function Connection = multiRecSummarize(varargin)
% CONNECTION = MULTIRECSUMMARIZE analyzes multi-recordings data and makes figures.
%  written by Taro Kiritani 10/15/2010 tarokiritani2008@u.northwestern.edu

switch varargin{1,1}
    case 'sp2sp'
        addpath('C:\Data\Taro\ANALYSIS\multi_rec_project\spinal2spinal')
        list = listSp2Sp;
    case 'st2sp'        
        addpath('C:\Data\Taro\ANALYSIS\multi_rec_project\striatal2spinal')
        list = list_St2Sp;        
    case 'st2st'
        addpath('C:\Data\Taro\ANALYSIS\multi_rec_project\striatal2striatal')
        list = listSt2St;        
end

summaryFigure = figure;
set(gco,'Name',varargin{1,1})
Connection = getEPSPdata(list);
% add cell type field in CONNECTION.
for k = 1:length(Connection)
    Connection(k).celltype = varargin{1,1};
end
showPostsynaptic(summaryFigure, Connection, list);
figtitle = title(varargin{1,1});
set(figtitle,'FontSize', 18)
Connection = subtractBaseline(summaryFigure, Connection, list);
Connection = normalize2firstEPSP(summaryFigure, Connection, list);
Connection = epspSummation(summaryFigure, Connection,list);
Connection = normEPSPplot(summaryFigure,Connection,list);
Connection = getBase2Troughs(summaryFigure,Connection,list);
[timeconstant,Connection] = plsticityAnalysis(Connection,summaryFigure,list);
Connection = measureIntegral(Connection);
% scatter plot of EPSP amplitude vs. ppr.
pprPlot(summaryFigure,Connection,list)

figure(summaryFigure);
subplot(449)
averageTrace(gca,summaryFigure,[4,4,13])

subplot(446)
averageTrace(gca,summaryFigure,[4,4,7],'on')

subplot(4,4,10)
averageTrace(gca,summaryFigure,[4,4,11],'on')
end

function Connection = getEPSPdata(list)
     for i = 1:length(list)
        hfig = open([list(i).filename]);
        lineObj = findobj(hfig,'Type','line');    
    
        % It is necessary to choose rigth traces. The order of line object does 
        % not correspond to the order of the lines.
        Connection(i).XData = get(lineObj(length(lineObj)+1-list(i).pre2post(1)),'XData');
        Connection(i).YDatapre = get(lineObj(length(lineObj)+1-list(i).pre2post(1)),'YData');
        Connection(i).YDatapost = get(lineObj(length(lineObj)+1-list(i).pre2post(2)),'YData');
        % make sure that the trace length is 1 sec (10000 points) long.
        Connection(i).XData = take1sec(Connection(i).XData);
        Connection(i).YDatapre = take1sec(Connection(i).YDatapre);
        Connection(i).YDatapost = take1sec(Connection(i).YDatapost);
        Connection(i).name = list(i).filename
        close
    end
end
function showPostsynaptic(summaryFigure, Connection, list)
figure(summaryFigure)
subplot(441)
for j = 1:length(Connection)
    hlines_raw(j) = plot(Connection(1).XData,Connection(j).YDatapost);
    set(hlines_raw(j),'Tag',list(j).filename)
    hold on
end
end
function Connection = subtractBaseline(summaryFigure, Connection, list)
figure(summaryFigure)    
subplot(445)
for k = 1:length(Connection)
    Connection(k).YDataBaselineSubtracted = Connection(k).YDatapost...
        -mean(Connection(k).YDatapost(1000:1999));
    hlines_baseline_subtracted(k) = plot(Connection(1).XData,Connection(k).YDataBaselineSubtracted);
    set(hlines_baseline_subtracted(k),'Tag',list(k).filename)
    hold on
end
title('baseline subtracted')
end
function Connection = normalize2firstEPSP(summaryFigure, Connection, list)
figure(summaryFigure)
subplot(449)
for k = 1:length(Connection)
    Connection(k).YDataNormalized = Connection(k).YDataBaselineSubtracted ./ max(Connection(k).YDataBaselineSubtracted(2000:2500));
    hlines_normalized(k) = plot(Connection(1).XData,Connection(k).YDataNormalized);
    Connection(k).snRatio = 1/std(Connection(k).YDataNormalized(1000:1999)) % signal is 1 because of the normalization.
    set(hlines_normalized(k),'Tag',list(k).filename)
    hold on
end



end
function Connection = epspSummation(summaryFigure, Connection,list)
% CONNECTION = EPSPSUMMATION(SUMMARYFIGURE,CONNECTION,LIST) measures the
% epsp amplitude of each epsp response from the baseline and then plot the
% change. It also adds YDATABASELINESUBTRACTED field to CONNECTION.

figure(summaryFigure)
subplot(442)
for m = 1:length(Connection)
    Connection(m).EPSP = [max(Connection(m).YDataBaselineSubtracted(2000:2499)),... 
        max(Connection(m).YDataBaselineSubtracted(2500:2999)),...
        max(Connection(m).YDataBaselineSubtracted(3000:3499)),...
        max(Connection(m).YDataBaselineSubtracted(3500:4000)),...
        max(Connection(m).YDataBaselineSubtracted(8500:8999))];
    Connection(m).firstEPSPamp = Connection(m).EPSP(1);
    hlines_EPSP(m) = plot([0 50 100 150 650],Connection(m).EPSP); 
    set(hlines_EPSP(m),'Tag',list(m).filename)
    hold on
end
end
function Connection = normEPSPplot(summaryFigure,Connection,list)
% NORMEPSPPLOT(SUMMARYFIGURE, CONNECTION, LIST) plots the summation of EPSP
% on SUMMARYFIGURE. EPSP amplitude is normalized to the first EPSP.
    figure(summaryFigure)
    subplot(446)
for n = 1:length(Connection)
    Connection(n).EPSP_normalized = [Connection(n).EPSP ./ Connection(n).EPSP(1)];
    hlines_EPSP_normalized(n) = plot([0 50 100 150 650],Connection(n).EPSP_normalized); 
    set(hlines_EPSP_normalized(n),'Tag',list(n).filename);
    hold on
end
end
function [timeconstant, Connection] = plsticityAnalysis(Connection,summaryFigure,list)
%     [TIMECONSTANT, CONNECTION] =
%     PLASTICITYANALYSIS(CONNECTION,SUMMARYFIGURE,LIST) calculates the
%     timeconstant for epsp decay. It also plots the change in epsp
%     amplitude. The decay is fitted with exponential. Then, the
%     exponential is subtracted from the data in order to measure the size
%     of each epsp.
timeconstant = zeros(length(Connection),4);
    
for p = 1:length(Connection)
    FittedRange = {[2200 2499],[2700 2999],[3200 3499]};
    Connection(p).EPSP_baseline_subtracted(1) = Connection(p).EPSP(1);
    Connection(p).EPSP_baseline_subtracted_norm(1) = 1;
    for q =1:3
        % find out the fitting exponential for the data.
        [A, lambda] = fitExpBound([1:300],...
            Connection(p).YDataBaselineSubtracted(FittedRange{1,q}(1):FittedRange{1,q}(2)));
        % the unit for A(2) is the [1/sec] (the reciprocal of XData)
        
        %f = @(x) A(1)*exp(-A(2)*x); % if two parameters are fitted.
        amp = Connection(p).YDataBaselineSubtracted(FittedRange{1,q}(2));
        % make sure that parameters are valid.
        if amp > 0 && A > 0
            f = @(x) amp*exp(-A * (x-300)); % if only one parameter (time constant) is fitted.
        else % linear regression is used if the parameters are invalid.
            param = polyfit([201:300],Connection(p).YDataBaselineSubtracted(FittedRange{1,q}(2) - 99:FittedRange{1,q}(2)),1);
            f = @(x) param(1) * x +param(2);
        end

        figure
        plot([1:300],Connection(p).YDataBaselineSubtracted(FittedRange{1,q}(1):FittedRange{1,q}(2)))
        hold on;
        ezplot(f,[1:300])
        timeconstant(p,q) = .1/A; % the unit now should be [ms]

        % subtract the decaying baseline.
%         decayFunc = @(x) A(1)*exp(-A(2)* (x - FittedRange{1,q}(1)/10000)); % this is the decalying baseline
         Connection(p).EPSP_baseline_subtracted(q+1) = ...
             max(Connection(p).YDataBaselineSubtracted(FittedRange{1,q}(2):FittedRange{1,q}(2)+500)...
         -f(1:501));
         Connection(p).EPSP_baseline_subtracted_norm(q+1) = ...
             Connection(p).EPSP_baseline_subtracted(q+1) / Connection(p).EPSP_baseline_subtracted(1);
         display(A)
 
%         figure;
%         plot(Connection(1).XData(2200:2999),Connection(p).YDataBaselineSubtracted(2200:2999)...
%         -decayFunc(Connection(1).XData(2200:2999)))
    end
    Connection(p).EPSP_baseline_subtracted_norm(5) = Connection(p).EPSP_normalized(5);
    
end

figure(summaryFigure);
subplot(4,4,10)
for r = 1:length(Connection)
     hEPSPsub = plot(Connection(r).EPSP_baseline_subtracted ./ Connection(r).EPSP_baseline_subtracted(1));
     set(hEPSPsub,'Tag',list(r).filename);
     hold on;
     
 end
end
function pprPlot(fig,Connection,list)
% PPRPLOT(FIG,CONNECTION) plot amplitude vs ppr plot on FIG. the 
    
for i = 1:length(Connection)
    ppr(i) = Connection(i).EPSP_baseline_subtracted(2)/Connection(i).EPSP_baseline_subtracted(1);
    amp(i) = Connection(i).EPSP_baseline_subtracted(1);
end

figure(fig)
subplot(4, 4, 14)

for j = 1:length(amp)
     hppr = scatter(amp(j),ppr(j));
     set(hppr,'Tag',list(j).filename)
     hold on
end

    xlabel('amplitude (\it{mV}\rm)')
    ylabel('PPR (2nd EPSP / 1st EPSP)')
end
function averageTrace(varargin)
% FUNCTION AVERAGEZTRACE(AXENUM, SUMMARYFIGURE, PLOTPOSITION) takes all the
% traces in AXENUM and then plot the average on SUMMARYFIGURE. The position
% of the axes is indicated by PLOTPOSITION.
% FUNCTION AVERAGEZTRACE(AXENUM, SUMMARYFIGURE, PLOTPOSITION, 'on') plots
% the error bar (sem) as well.

axesnum = varargin{1,1};
summaryFigure = varargin{1,2};
plotPosition = varargin{1,3};

traceHandles = get(axesnum,'Children');
Trace = get(traceHandles,'YData');
TraceMean = mean(cell2mat(Trace));
figure(summaryFigure)

if isvector(plotPosition)
    subplot(plotPosition(1),plotPosition(2),plotPosition(3))
else
    subplot(plotPosition)
end

if nargin >3
    flag = varargin{1,4};
    if flag == 'on'
       TraceSem = std(cell2mat(Trace)) / sqrt(length(traceHandles));
       errorbar(TraceMean,TraceSem) 
    end
else
    plot(TraceMean)
end
end
function trimmedTrace = take1sec(trace)   
    if length(trace) > 10000
        trimmedTrace = trace(1:10000);
    else
        trimmedTrace = trace;
    end
end
function Connection = getBase2Troughs(summaryFigure,Connection,list)
% FUNCTION CONNECTION = GETBASE2TROUGHS(SUMMARYFIGURE,CONNECTION,LIST)
% measures the amplitude of troughs from the baseline.
for k = 1:length(Connection)
    Connection(k).TroughsAmp = Connection(k).YDataNormalized([2499,2999,3499]);
end
end
function Connection = measureIntegral(Connection)
% FUNCTION CONNECTION = MEASUREINTEGRAL(SUMMARYFIGURE,CONNECTION,LIST)
% measures the area under curve of EPSP.
integralWindow = 40; % the unit is ms.

% start and end points for integration should take into account the action
% potential timing.
for k = 1:length(Connection)
    aptime = Connection(k).YDatapre;
    [ampmax aptime] = max(aptime(1:2500));
    Integral = Connection(k).YDataBaselineSubtracted(aptime:aptime + integralWindow * 10); 
    Connection(k).Integral = sum(Integral) / 10 ; % unit is (mV * ms)
end

% when area under curve is measured after the last AP, baseline should be
% shifted, to cmpensate the deviation.
for kk = 1:length(Connection)
    aptime = Connection(kk).YDatapre;
    [ampmax aptime] = max(aptime(8500:9000));
    aptime = aptime + 8499;
    Integral = Connection(kk).YDataBaselineSubtracted(aptime:aptime + integralWindow * 10)...
        - mean(Connection(kk).YDataBaselineSubtracted(8400:8500));
    Connection(kk).Integrallate = sum(Integral) / 10 ; % unit is (mV * ms)
end
end
